Low temperature thermal conductivity in a d-wave superconductor with coexisting 
charge order: Effect of self-consistent disorder and vertex corrections 



00 
O 

o 

(N 

> 
O 



o 
o 

u 

Oh: 



I 

O 

o 



> 
o 



oo 
O 



X 



Philip R. Schiff and Adam C. Durst 
Department of Physics and Astronomy, Stony Brook University, Stony Brook, NY 1 1794-3800, USj^ 

(Dated: November 1, 2008) 

Given the experimental evidence of charge order in the underdoped cuprate superconductors, we 
consider the effect of coexisting charge order on low-temperature thermal transport in a d-wave 
superconductor. Using a phenomenological Hamiltonian that describes a two-dimensional system in 
the presence of a Q = (tt, 0) charge density wave and d-wave superconducting order, and including 
the effects of weak impurity scattering, we compute the self-energy of the quasiparticles within the 
self-consistent Born approximation, and calculate the zero-temperature thermal conductivity using 
linear response formalism. We find that vertex corrections within the ladder approximation do not 
significantly modify the bare-bubble result that was previously calculated. However, self-consistent 
treatment of the disorder does modify the charge-order-dependence of the thermal conductivity 
tensor, in that the magnitude of charge order required for the system to become effectively gapped 
is renormalized, generally to a smaller value. 

PACS numbers: 74.72-h, 74.25.Fy 
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I. INTRODUCTION 

The superconducting phase of the cuprate supercon- 
ductors exhibits d-wave pairing symmetry.^ As such, 
there exist four nodal points on the two-dimensional 
Fermi surface at which the quasiparticle excitations are 
gapless, and quasiparticles excited in the vicinity of a 
node behave like massless Dirac fermions.— The 
presence of impurities enhances the density of states 
at low energy? resulting in a universal limit (T — > 
0, f2 — > 0) where the thermal conductivity is indepen- 
dent of disorder i^'^'^'^ d°did^ Calculations have shown 
that the thermal conductivity retains this universal char- 
acter even upon the inclusion of vertex corrections.^^ Ex- 
periments have confirmed the validity of this quasipar- 
ticle picture of transport by observing their universal- 
limit contribution to the thermal conductivity, and 
thereby measuring the anisotropy of the the Dirac nodes, 
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For some time, there has been significant 
interes1i2ii2^iS£i22i2^ in the idea of additional types 
of order coexisting with c?-wave superconductivity 
(dSC) in the cuprates. And in recent years, as the 
underdoped regime of the phase diagram has been 
explored in greater detail, evidence of coexisting order 
has grown substantiallyi2i. Particularly intriguing 
has been the evidence of checkerboard charge order 
revealed via scanning tunnelling microscopy (STM) 
experiment: . ^^•^"•'^^•^^•'^^•^'^•^^•^^•^'^■^^•^^•"'^"•'^"'^•'^^ 

And if charge order coexists with d-wave superconduc- 
tivity in the underdoped cuprates, it begs the question 
of how the quasiparticle excitation spectrum is modified. 
Previous work^ has shown that even with the addition 
of a charge or spin density wave to the dSC hamiltonian, 
the low-energy excitation spectrum remains gapless as 
long as a harmonic of the ordering vector does not nest 
the nodal points of the combined hamiltonian. However, 



if the coexisting order is strong enough, the nodal points 
can move to fc-space locations where they are nested by 
the ordering vector, at which point the excitation spec- 
trum becomes fully gapped^ii^i^ 

Such a nodal transition should have dramatic conse- 
quences for low-temperature thermal transport, the de- 
tails of which were studied in Ref. \^ That paper con- 
sidered the case of a conventional s-wave charge den- 
sity wave (CDW) of wave vector Q = (tt, 0) coexist- 
ing with c?-wave superconductivity. It showed that the 
zero-temperature thermal conductivity vanishes, as ex- 
pected, once charge order is of sufficient magnitude to 
gap the quasiparticle spectrum. In addition, the depen- 
dence of zero-temperature thermal transport was calcu- 
lated and revealed to be disorder-dependent. Hence, in 
the presence of charge order, the universal-limit is no 
longer universal. This result is in line with the results of 
recent measurement a"'^^i"^^'^°i^^'^^i^^'^"^ of the underdoped 
cuprates, as well as other calculation a^^'^^ . 

We extend the work of Ref. 113 herein. We consider the 
same physical system, but employ a more sophisticated 
model of disorder that includes the effects of impurity 
scattering within the self-consistent Born approximation. 
We find that this self-consistent model of disorder re- 
quires that off-diagonal components be retained in our 
matrix self-energy. These additional components lead to 
a renormalization of the critical value of charge order 
beyond which the thermal conductivity vanishes. Fur- 
thermore, we include the contribution of vertex correc- 
tions within our diagrammatic thermal transport calcu- 
lation. While vertex corrections become more important 
as charge order increases, especially for long-ranged im- 
purity potentials, we find that for reasonable parameter 
values, they do not significantly modify the bare-bubble 
result. 

In Sec. [ni we introduce the model hamiltonian of the 
dSC-|-CDW system, describe the effect charge ordering 
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has on the nodal excitations, and present our model for 
disorder. In Sec. IIII Ai a numerical procedure for com- 
puting the self-energy within the self-consistent Born ap- 
proximation is outlined. The results of its application 
in the relevant region of parameter space are presented 
in Sec. IIII Bl In Sec. llVi we calculate the thermal con- 
ductivity using a diagrammatic Kubo formula approach, 
including vertex corrections within the ladder approxi- 
mation. An analysis of the vertex-corrected results and a 
calculation of the clean-limit thermal conductivity is pre- 
sented in Sec.|Vl Also in this section, we discuss how our 
self-consistent model of disorder renormalizes the nodal 
transition point, the value of charge order parameter at 
which the nodes effectively vanish. Conclusions are pre- 
sented in Sec. IVll 



II. MODEL 

We employ the phenomenological hamiltonian of 
Ref. 113 in order to calculate the low-temperature thermal 
conductivity of the fcrmionic excitations of a d-wave su- 
perconductor with a Q = (tt, 0) charge density wave, in 
the presence of a small but nonzero density of point-like 
impurity scatterers. The presence of d-wave supercon- 
ducting order contributes a term to the hamiltonian 



Hdsc = t: 



ka 



h.c. 



(1) 



where is a typical tight-binding dispersion, and 
an order parameter of d^2_y2 symmetry. Due to the d- 
wave nature of the gap, nodal excitations exist in the 
(±7r, ±7r) directions with respect to the origin. The loca- 
tions of these nodes in the absence of charge ordering are 
close to the points (±7r/2, ±7r/2), and are denoted with 
white dots in Fig. [1] These low energy excitations are 
massless anisotropic Dirac fermions. That is, the elec- 
tron dispersion and pair function are linear functions of 
momentum in the vicinity of these nodal locations. We 
will refer to the slopes of the electron dispersion and pair 
function, defined by v/ = ^ and va = as the 

Fermi velocity and gap velocity respectively. The en- 
ergy of the quasiparticles in the vicinity of the nodes is 

given by Ek = \Jv'jk'l -f u^fc|, where ki and /c2 are the 

momentum displacements (from the nodes) in directions 
perpendicular to and parallel to the Fermi surface. The 
universal-limit (T — s- 0, f2 — s- 0) transport^properties of 
these quasiparticles was explored in Ref. il2i. 

While experiments have revealed evidence of a number 
of varieties of spin and charge order, the system described 
in this paper will be restricted to the addition of a site- 
centered charge density wave of wave vector Q = (tt, 0), 
which contributes a term to the hamiltonian 



HcDW = '^a.kc\^Ck+Qa + h.c. 



(2) 



ka 



The charge density wave doubles the unit cell, reduc- 
ing the Brillouin zone to the shaded portion seen in 
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FIG. 1: Illustrated is the Brillouin zone for our model, reduced 
to the shaded region by unit-cell-doubling charge order. The 
^ = nodal locations are illustrated by white dots. They are 
displaced by a distance fco from the (±-|,±-|) points (stars). 
As the charge density wave's amplitude increases, the location 
of the gapless excitations evolves along curved paths toward 
the (±-|, ±-|) points, until i/) reaches ipc., when the spectrum 
becomes gapped because the nodes are nested by the charge 
density wave-vector. The gray dots depict the images of the 
nodes in the second reduced Brillouin zone. 



Fig. [TJ Restricting summations over momentum space 
to the reduced Brillouin zone, and invoking the charge 
density wave's time-reversal symmetry and commensu- 
rability with the reciprocal lattice, we are able to write 
the hamiltonian as 



dSC 
k 



CDW 



(3) 



where 



Hk = 



Afe 

Vo 



Afe 

-Cfe 






(4) 



Ak+Q 

is a matrix in the basis of extended-Nambu vectors, 



/ Ck, \ 



CfeT 
t 

Cfe-I-Qt 
V'^-fe-QJ, 



(4t 



-kl Cfc+Qi C-k-Ql 



(5) 

and ip represents the constant value taken at the nodes by 



^fe+Q- 



the charge density wave order parameter — ak - 

The onset of the charge order modifies the energy spec- 
trum of the clean hamiltonian so that the locations of the 
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nodes evolve along curved paths towards the (±^,±5) 
points at the edges of the reduced Brillouin zone, as was 
noted in Ref. |4J. "Ghost" nodes, their images in what is 
now the second reduced Brillouin zone, evolve the same 
way, until the charge density wave is strong enough that 
the nodes and ghost nodes collide at those (±7r/2, ±7r/2) 
points. When that occurs, Q nests two of the nodes, gap- 
ping the spectrum so that low temperature quasiparticle 
transport is no longer possible. We define the value of 
at which this occurs as ipc- Due to the nodal properties 
of the quasiparticles, all functions of momentum space k 
can be parametrized in terms of a node index j, and local 
coordinates pi and p2 in the vicinity of each node. We 
choose to parametrize our functions using symmetrized 
coordinates centered at (±7r/2, ±7r/2), 

Efc = i'c + Ppi Afe = 

P 

Cfc+Q V^c + I3p2 ^k+Q = -j^Pi (6) 

where we have rescaled y/vJvXki = pi for the coordi- 
nate normal to Fermi surface, ^w/ua ^2 = P2 for the 
coordinate parallel to Fermi surface, and introduced the 

definition /3 = \J'^- In this coordinate system, the dis- 
placement of the original node locations from the collision 
points is given hy ill c A sum over momentum space is 
therefore performed by summing over nodes, and inte- 
grating over each node's contribution, as follows. 

1 4 rPa rPa 

= E/ dpi/ dp2 f'^^HPl,P2) (7) 

where the factor of ^ comes from extending the integrals 
to all pi and p2 , rather than just the shaded part depicted 
in Fig. [Tl and po is a high-energy cutoff. 

At sufficiently low temperatures, the thermal con- 
ductivity is dominated by the nodal excitations, since 
phonon modes are frozen out, and other quasiparticles 
are exponentially rare. Using this fact, we can calculate 
the low temperature thermal conductivity of the system 
using linear response formalism. 

We incorporate disorder into the model by including 
scattering events from randomly distributed impurities. 
Because the quasiparticles are nodal, only limited infor- 
mation about the scattering potential is needed, in par- 
ticular, the amplitudes Vi, V2 and V3, for intra-node, ad- 
jacent node, and opposite node scattering respectively, 
as explained in Ref. [ij. We calculate the thermal con- 
ductivity using linear response formalism, wherein we 
obtain the retarded current-current correlation function 

by analytic continuation of the corresponding Matsubara 
correlator.5Zi58 ^ 

In Ref. 113, using a simplified model for disorder, where 
the self-energy was assumed to be a negative imaginary 
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FIG. 2: Feynman diagram depicting self-energy in the self- 
consistent Born approximation. The double line represents 
the dressed propagator, the dashed line represents the interac- 
tion with the impurity, and the cross represents the impurity 
density. 

scalar, the thermal conductivity was calculated as a func- 
tion of ip, and found to vanish for ip > ipc- We now 
improve upon that result by calculating the self-energy 
within the self-consistent Born approximation, and by 
including vertex corrections within the ladder approxi- 
mation in our calculation of the thermal conductivity. 



III. SELF-ENERGY 

A. SCBA Calculation 

Within the self-consistent Born approximation 
(SCBA), the self-energy tensor is given by 

E(k,w) ^ nin-,p^\Vkk'f {(To <^ T3)g{k,uj){aQ 1^ T3) (8) 

k' 

where riimp is the impurity density and Vkk' = 
Vfcfc'(cro (8) T3) accompanies each scattering event, as seen 
in Fig. [21 The tilde signifies an operator in the extended- 
Nambu basis, and the cr's and t's are Pauli matrices 
in charge-order-coupled and particle-hole spaces respec- 
tively. t/(k, oj) is the full Green's function, whose relation 
to the bare Green's function Qo(k,oj) and the self-energy 
S(k, is given by Dyson's equation 

g{Ku;)^{g^\k,u;)-^{k,u;))-\ (9) 

the bare Green's function having been determined by 

go{k,u) = {ujl-Hk)-\ (10) 

Eq. ([8]) and Eq. ^ define a set of integral equations 
for the self-energy E(k, lj). For the calculation of the 
universal-limit thermal conductivity, it is sufficient to 
find the zero-frequency limit of the self-energy. In its 
present form, E has 32 real components. Below, we 
demonstrate that this number can be reduced further to 
six components. 
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If we write the Green's function as 
1 I Ga Gb 



^den 



Qc Qd I ' 



where 



(11) 



(12) 



i=0 



then the self-energy can be written as the set of 16 com- 
plex equations (for a ~ {A, B, C, D}, i — {0, 1, 2, 3}) 



> ^2 1^ '^T^Ga 



if den 



2 GaiiPl,P2) 



Gden (Pl,P2) 



(13) 



where 



ni{V{'+2Vi + Vi) 



and the 



+1, i = 0,3 
-1, z-1,2 

final line is realized by using the notation of Eq. ([6]) and 
Eq. ([7]) and completing the sum over nodes. From the 
symmetries of the hamiltonian, we are able to ascertain 
certain symmetries the bare Green's function will obey, 
specifically. 



(14) 



GAoiP2,Pl) 




{Pl,P2) 


GAliP2,Pl) 


- c(") 


{Pl,P2) 


GAiiP2,Pl) 


- 


{Pl,P2) 


^50(^2, Pi) 


- 


{Pl,P2) 


^^Bl(P2,Pl) 


- 

— y^i 


{Pl,P2) 


^B2(P2,Pl) 


- c(o) 

— yc2 


{Pl,P2) 


^B3(P2,Pl) 


- c(o) 


{Pl,P2) 


GdliP2,Pl) 


= 

■^clcn 


{Pl,P2) 





— ^Di 






= ^Ci 


i = 0,1,2,3 


Sb2 


= SC2=0 





In addition, the realization that the integration is also 
symmetric with respect to exchange of pi andp2, coupled 
with these symmetries, lead to relations for self-energy 
components 



(15) 



so that we see a reduction from 32 components of 
the self-energy to 6 independent components:{I]Qi} = 
{S^o, Syii, Syia, Ebo, Sbi, Ssa}- A self-consistent self- 
energy must therefore satisfy 6 coupled integral equations 
given by Eq. p^ . 

The self-consistent calculation of the self-energy pro- 
ceeds by applying the following scheme: First, a guess 
is made as to which self-energy components will be in- 
cluded. The full Green's function corresponding to such 
a self-energy is then obtained from Dyson's equation, 
Eq. The quantitative values of the E^i's are then 



determined as follows: An initial guess for the quanti- 
tative values of each of the S^i's is made, and the six 
integrals of Eq. (|13p are computed numerically, which 
provides the next set of guesses for {Eai}. This process 
is repeated until a stable solution is reached. Finally, the 
resulting solutions must be checked that they are consis- 
tent with the initial guess for the form of E. If they are, 
the self-consistent calculation is complete. 

We begin with the simplest assumption, that — 

— iro((To tq), where Fq is the zero- frequency limit of the 
scattering rate. The superscript indicates that this is the 
first guess for E. The Green's function components are 
computed, which gives the explicit form of Eq. ([8]). Upon 
evaluating the numerics, it is seen that this first iteration 
generates a nonzero (real and negative) term for Esi. So, 
the diagonal self-energy assumption turns out to be in- 
consistent, in contrast to the situation for ip ~ 0. We 
then modify our guess, assuming self-energy of the form 

^ -iFo(cro ® To) - Bi{(Ti (g) Ti). The Green's func- 
tion is computed again, using Dyson's equation, and the 
self-energy equations are obtained explicitly. It is noted 
that the symmetries of Eq. ([H]) still hold. Again, the 
equations (fT3)) are solved iteratively; the result is a non- 
zero Es3 component as well. Once again, the Green's 
functions are modified to incorporate this term, and the 
iterative scheme is applied. Calculation of the self-energy 
based on the assumption 



E(3) = -iVoia^To) - Bii<^ 



Ti) ~ Bsiai (g) T3) 

ro,Bi,S3 >o(i6) 



generates Fo, Bi, and B3 that are much larger than any 
remaining terms, and hence provides the self-consistent 
values of E^O; E^i and Ess- A plot of the 6 components 
of E is displayed in Fig. [3] for a representative parameter 
set, where we see that the three terms of the ansatz are 
indeed dominant. For the remainder of this paper, the 
effect of the E^n, T,a3 and Esq components will be ig- 
nored. The self-consistent Green's functions are provided 
in Appendix [XI while additional details of the self-energy 
calculation are discussed in Appendix [B] 



B. SCBA Results 

In order to discuss the numerical results contained in 
this paper, it is necessary to make a note about the units 
employed. The following discussion of units applies as 
well to the numerical analysis of the results of the thermal 
conductivity calculation in Sec.|Vl Because we are study- 
ing the evolution of the system with respect to increasing 
CDW order parameter ip, we wish to express energies in 
units of the value of which gaps the clean system. In 
order to do this, the cutoff po is fixed such that the Bril- 
louin zone being integrated over in Eq. ([7]) has the correct 
area. In this way, pq sets the scale of the product VfVA, 

a parameter /3 



is defined to represent the veloc- 
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FIG. 3: Components of self-energy computed using iter- 
ative procedure described in Sec. Hira The third itera- 
tion self-energy, S'^-*, is shown here. The dominance of 
Fo = -Im(EAo),Bi = -Re(EBi), and B3 = -Re(Es3) over 
other components establishes this third iteration as yielding 
the (approximately) self-consistent value of the self-energy. 
Eai and Ea3 overlap. 



ity anisotropy. Then, ^ — j^^/vJva, so that we may 
eliminate the frequently occurring parameter ATTVfV/^ by 
expressing lengths in units of -^a w 2.25a. Impurity 
density ni,„p is thus recast in terms of impurity frac- 
tions z according to njmp — —z. Finally, the parameters 
of the scattering potential are recast in terms of their 
anisotropy. We define V2 = R2V1 and V3 = R3V1. 

With these modifications, the original set of pa- 
rameters, Vi, 1^2, V3, w/, wajPoj "01 "^c} is reduced to 
{z, Vi, i?2i /SjPoj V'}. For the work contained herein, 
the cutoff po is fixed at po = 100. The self-energy in 
the self-consistent Born approximation was computed for 
different scattering potentials as a function of impurity 
fraction and CDW order parameter ■0. Since it was found 
that three of the components, Saoj ^bi and Ss3, dom- 
inate over the others, we will subsequently analyze only 
those three components, referring to their magnitudes as 
Fq, -Bi, and -B3 respectively. 

As z — > 0, the Green's functions become impossibly 
peaked from a numerical point of view. For sufhciently 
large z, depending on the strength of the scatterers, the 
Born approximation breaks down. Given a scattering 
strength of Vi = 110, cutoff po = 100, scattering poten- 
tials that fall off slowly in fc-space and velocity anisotropy 
ratios (3 = ^/vJJva — 1, 2, 3, 4, this puts the range of z in 
which our numerics may be applied at roughly between 
one half and one percent. 

Some results for £(7/;) , for several values of z, are shown 
in Figs. [4] and [5] These plots correspond to the same 
parameters, except that Fig. [J] illustrates the Vf — v\ 
case, and Fig. [5] illustrates Vf — 16wa. In all cases it is 
seen that 

Bi{^,z)~bi{z)il; 

B3{^,z)c^b3{z)^ (17) 



0.15 F 




FIG. 4: Effect of disorder on charge-order-dependence of self- 
energy components. To satisfy Dyson's equation, it is neces- 
sary to include three (extended-Nambu space) components of 
the self energy. Their self-consistent values are plotted here 
for several different values of impurity fraction z. Here, the 
scattering potential is given in our three parameter model as 
{Vi, R2, R3} = {110, 0.9, 0.8}, which represents a fairly short- 
ranged potential. These results are for the case of isotropic 
nodes {vf = va)- All energies are in units of t/jc. 

where the dependence of B3, bi, and 63 on the re- 
maining parameters is implicit. For much of the parame- 
ter space sampled, Fq does not have much ip dependence, 
except that it typically rises and then falls to zero at some 
sufficiently large ip < This feature will be revisited in 
Sec. |Vl wherein it is explained that this vanishing scat- 
tering rate coincides with vanishing thermal conductivity, 
and corresponds to the point at which the system be- 
comes effectively gapped and our nodal approximations 
break down. The value of ip at which this occurs depends 
on the entire set of parameters used, and will be referred 
to as Ip*. The observed z dependence is not very sur- 
prising, in light of Eq. (|13p. The self-energy components 
depend on z roughly according to 

Fo poexp(-i) 
z 

El ^ z 

B3 ^ z (18) 

as can be seen in Fig. [6l When ip = 0, Tq i s gi ven by 
the closed-form expression obtained in Ref. Il2l . Fq = 

Poexp(:^), where c = "'^^o t?,^^ """^^ ^ . For finite ih, 



6 




FIG. 5: Effect of disorder on charge-order-dependence of self- 
energy components. This figure illustrates the case where 
Vf = l&VA- The scattering potential is again given by 
{Vi,R2,Ri,} ~ {110,0.9,0.8}, representing a fairly short- 
ranged potential. The plots for B\ and Bs terminate before 
ijj reaches V'c because for sufficiently large V'l the excitations 
become gapped and our nodal approximations break down. 
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FIG. 6: Effect of charge order on disorder-dependence of self- 
energy components. Nonzero components of are shown 
for impurity fraction z ranging from 0.5 to 1.0%, for charge 
order parameter ^/j^O, 0.2, 0.4, 0.6, and 0.8 (in units of i/)c). 
These results are for scattering parameters {Vi,i?2,J?3} = 
{110,0.9,0.8} and «/ — va- Similar results are obtained for 
the case of anisotropic nodes. 



this precise form does not hold, but the strong z de- 
pendence of To remains, in contrast to that of Bi and 
B3. Note that the z dependence of Bi and B3 is roughly 
linear for <C V'*- As V' approaches V'* the functions 
diverge slightly from linearity. Results for several values 
oi ip < ip* are shown in the figure. 



where a generalized velocity is defined as 

= {as (E) T„) My = (ao ® r„) 



(21) 



IV. THERMAL CONDUCTIVITY 



Thermal conductivity was calculated using the Kubo 
formula^i^. 



K{n,T) _ ImnRet(») 



(19) 



where nRet(f^) is the retarded thermal current-current 
correlation function. To find this correlator, it is nec- 
essary to first compute the appropriate thermal current 
operator. For our model hamiltonian, this is done in Ref. 
I47I with the result 



(v/M + vam) tpk+q, (20) 



where a — {/, A} and — {t3,ti} for Fermi and gap 
velocities respectively. 

To calculate a thermal conductivity that satisfies Ward 
identities, vertex corrections must be included on the 
same footing as the self-energy corrections to the single 
particle Green's function. The details of this calculation 
are similar to those performed in Appendix B of Ref. [T^ . 
The impurity scattering diagrams which contribute to 
the ladder series of diagrams are included by express- 
ing the correlation function in terms of a dressed vertex, 
as shown in Fig. 7(a) The current-current correlation 



function is obtained from this dressed bubble. The bare 
current operator of Eq. (jSO]) is associated with one ver- 
tex of the bubble, while the dressed vertex of Fig. |7(b) 
is associated with the other. Evaluating Fig. 7(a)[ we 
find that the current-current correlation function takes 
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FIG. 7: (a) Feynman diagram representing the correlation 
function H™ in terms of a bare vertex j™, and a dressed 
vertex FJ^. (b) Feynman diagram representing the (ladder 
series) dressed vertex in terms of the bare vertex and the 
Born scattering event. 



the form 



a,/3=/,A 



Tr 



GiVcK'M^g^vpM-T^J (22) 



where Qi = t/(k, zw), G2 = ^(k, + ifi), and = 
r^(k, iw, ifi) represents the dressed vertex depicted in 
Fig. |7(b)[ The Greek indices denote "Fermi" and "gap" 
terms, while the Roman indices denote the position space 
components of the tensor. We use Fig. |7(b)| to find the 
form of the vertex equation, and then make the ansatz 
that 



(23) 



rf3{k,iuj,in) = {t + A{\k.\,iLu,in)]k, 
which leads to the scalar equations 

f'^{k,iuj,in) = fc„(i + a;,'). 



(24) 



Looking for solutions of this form, we see that the scalar 
vertex function is 

^ ^ ^ ^ ^ ^ )^'" 

= n,J2M^'Vkk'g2M^{t + A''^)giVk'k-^. (25) 



■/3 



Since we are working with nodal quasiparticles, we uti- 
lize the parametrization of Eq. ([7]), so that the vertex 



function is now a function of node index j and local 
mentum p 



mo- 



A^ = n; 



—JJ —J 0^ ,(]) ' J ST^2y 



M^{a7^3)G2MJ^\l + A'/,)gi{<j7^3)- (26) 
Arbitrarily choosing j = 1, then for j' = {1, 2, 3, 4} 
L(i') yfcO") 

= n -1 -1 11 = 111-1 -11 

= {1,-1, -1,1} ^ = {1,1,-1,-1}. (27) 

'^2x ^2y 

Using the node space matrix representing the 3- 
parameter scattering potential 



V 

—33 



(Vx V2 V2^ 

V2 Vi V2 V3 

V3 V2 Vi V2 

\V2 V3 V2 Vi, 



(28) 



we obtain for the vertex equation 

j ^M^{<^T3)g2M^(i + A^)ei(a^3) 



(29) 



^2 ^2 

where 7 = n\mr>-^ —■ The correlator then becomes 

iuj k 

Tr(^iM^^2M^(l+A^)) 
= l^(z..+ !^^)2^(fc«fcg) 

iui j=l 

d^p 



Since 



4 

E = 2 ((1 - 5o.p)rj„, + 5o.p) 5mn (31) 



we can write 



where 



mn _ 



and 



(32) 

(1 " 6af3)rim + SaP^Smn (33) 

/^2 
-^giM^g2Mi (34) 



1 VgVfl 
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To calculate the conductivity, we will need Tr(/^) and 
Tr(/^Ap. For if; — 0, it is possible to compute the 
integral in Eq. (|34p analytically, but for general ijj we 
had to compute the integrals numerically. We note that 
if we write 



I A Ib 

Ic Id 



(35) 



apply the symmetry properties of Eq. (I14|) and reverse 
the order of integration of pi and p2, then I a = Id, and 
Ib ^ Ic, so that the most general expansion of in 
Nambu space is 



Then 



Tr(/, 



mn \ 
af3 ) 



\ i=0 i'=0 

4(^q']3")oo, 



while if we use the same expansion for 

1 3 



i=0 i'=0 



we find 



Tr(C/A^) = E E (C/r)..'(A^) 
ij=o i'j'=a 



33 



Tr(o-iO-j 



= 4EE(CA^'(^^)^ 



Then Eq. ^ becomes 



4(A^),,- = Tr (^(a, ® T,OA^ 

d2„ , ~ , 



7Tr(L^,,,( 



1 + A'i' 



where 



L 



(3ii' 



IT 



(36) 



(37) 



(38) 



(39) 



(40) 
(41) 



(42) 



The symmetries of Q which were used to see which com- 
ponents of /™ were can also be applied to i^--, 
with the result that (L^,,,)a = (^^,,0^, (i^.^Os = 



ViiL0u')c, where = 



Since all that 



+1, i = 0,l 

^f^U')C, wn«r« = ^ _i^,^2,3 

is required for the conductivity is i = 0, 1, we use the 
expansion 



= E Y.(^'^3')iL0^^')3 

3=0j'=0 



(43) 



so that 



(A^).., = i7Tr(i^,,,(l + A^) 
= i7Tr(^^(i^,,,).y('^7^.') 

3=0j'=0 
1 3 

+ E E (-^/3"')jj'(Aj3)fefc') 
jfe=Oi'fe'=0 

1 3 
]=0j'=0 



(44) 



The thermal conductivity is obtained from the retarded 
current-current correlation function 



K™"(f7) _ 1 Im(n™(f7)) 



T 



T 



(45) 



where nict(f^) = n(zri ^ -I- i6). To get the retarded 
correlator we first perform the Matsubara summation. 
Consider the summand of Eq. 1321 which we redefine ac- 
cording to 



(46) 



J{iuj, iLU + ifl) = Tr + A^) 



The function J{iuj, iu+ifl) is of the form J{iuj, iuj+i^) =^ 
f {A{iuj) B {iu! + ift)) where A and B are dressed Green's 
functions of a complex variable z = iujn, so that J is 
analytic with branch cuts occurring where z and z + i^l 
are real. The Matsubara summation needed is performed 
by integrating on a circular path of infinite radius, so that 
the only contribution is from just above and just below 
the branch cuts, 



1 



n™" = -c,"- i dz(z + -)V(z, z + in) 



de7i/(e) 



(e + Y)^(>^(e + + i^) - J(e ~'i6,e + i^)) 
+ (e - Y)^(^(e + iS) - J(e -in,e- z<5))) . (47) 

To obtain the retarded function, we analytically continue 
iQ ^ fl + i6. Then we let e — > e + in the third and 
fourth terms, so that 



;Re( J^^«(e, e + 17) - J^f (e, 6 + f7))(48) 
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where J^^ and J^^ are defined by Eqs. (ge]) and ((44)) 
and are composed of the universal-hmit Green's functions 
given in Appendix [A] Taking the imaginary part, we find 



K'""(f7,T) 

T 



T 



l^C^ Re (J^^{e,e + n))-J^fie,e + n)). (49) 

a/3 

In taking the — * limit, the difference in Fermi 
functions becomes a derivative. Evaluating the integral. 



, we find that 



C"Re ( J^^«(0, 0) ~ J^f (0, 0)) .(50) 



T 

That K'^y = ^ is seen from Eq. Finally, 
since the a (3 integrals are traceless, the result for the 
thermal conductivity is 



T 



3 VfVA 8 



;^a;l'(0,0)- J5f(0,0)). (51) 



V. RESULTS 

For a discussion of the units employed in the analy- 
sis, one can refer to Sec. IIII Bl The reduced set of pa- 
rameters for the model is {z,Vi, R2, R3, P,po,i^}- We 
explored a limited region of this parameter space, cal- 
culating the integrals and solving the matrix equation 
numerically. In particular, we looked at the tp depen- 
dence of K. To vary the anisotropy of the scattering po- 
tential, we considered the {i?2,-R3} values of {0.9,0.8}, 
{0.7,0.6}, and {0.5,0.3}, and kept fixed the constant c 
(given after Eq. (I13p )bv appropriately modifying Vi. For 
{i?2,i?3} = {0.9,0.8}, we used Vi = 110. The rationale 
for keeping c fixed is that the self-energy depends only on 
c, P and pq. Additionally, we explored the dependence 
of the thermal conductivity on impurity fraction z and 
velocity anisotropy /3. For all computations we set the 
cutoff po — 100; this simply fixes a particular value of the 
product VfVA for these calculations. 
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FIG. 



0.4 , 0.6 



Vertex-corrected thermal conductivity, in units of 



the universal conductivity kq/T = ^{vf/v^ +VA/vf)- This 
data reflects a short range scattering potential {Vi, 7?2, JSs} = 
{110,0.9,0.8}, impurity fraction z—0.01, and isotropic Dirac 
quasiparticles {vf = «a). The inset displays the discrepancy 
between the bare-bubble and vertex-corrected results, in units 
of the bare-bubble result. It is clear that the vertex correc- 
tions are of little quantitative importance for these particular 
parameters. 




0.4 , 0.6 



1.0 



The importance of including the vertex corrections 
is determined by comparing the vertex corrected ther- 
mal conductivity with that of the bare-bubble. If 
— << 1 tor a region of parameter space, then m 
that regime the bare-bubble results can be used instead. 
This is of threefold practicality: the bare-bubble results 
are less computationally expensive, the bare-bubble ex- 
pression is much simpler to analyze, and other hamilto- 
nians could be more easily studied. 

The bare bubble thermal conductivity can be obtained 
by setting Ao ^ in Eq. or by using a spectral 



FIG. 9: Vertex-corrected thermal conductivity, in units of the 
universal conductivity kq/T = ■|f-(w//i'A +VA/vf). This fig- 
ure portrays the effect that a different scattering potential has 
on the importance of vertex corrections. Here, a longer range 
potential {Vi,-R2,-R3} ~ {140,0.5,0.3} was used, again with 
impurity fraction 2=0.01 and vf — va- The inset displays the 
discrepancy between the bare-bubble and vertex-corrected re- 
sults, in units of the bare-bubble result. From this, we deter- 
mine that vertex corrections make a more substantial correc- 
tion as the forward scattering limit is approached, but only 
once the charge ordering is quite strong. 
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FIG. 



10; Vertex-corrected thermal conductivity, 



m units 



of the universal conductivity kq/T = -^(vf/vA + v^/vf). 
Again, a short-ranged scattering potential, {Vi, R2, Ra} = 
{110,0.9,0.8} and isotropic nodes {vf = wa) are used. This 
figure displays the efi'ect of a smaller impurity fraction than 
that depicted in Fig|S] The inset displays the discrepancy be- 
tween the bare-bubble and vertex-corrected results, in units 
of the bare-bubble result; since the scattering potential falls 
off slowly (in k-space) here, the vertex corrections are again 
quite unimportant. 




0.4 0.6 



1.0 



FIG. 11: Vertex-corrected thermal conductivity, in units 
of the universal conductivity ko/T = ^{vf/vA + 
VA/vf), for short-ranged scattering potential, {Vi,-R2,-R3} ~ 
{110,0.9,0.8} and impurity fraction z — 0.01. These calcu- 
lations differ from those of Fig[8] in that they apply to the 
case of a more anisotropic Dirac spectrum with Vf = 9va- 
The thermal conductivity has a qualitatively similar t/j de- 
pendence, but vanishes for a smaller value of tp than for the 
isotropic case. The inset displays the discrepancy between 
the bare-bubble and vertex-corrected results, in units of the 
bare-bubble result; again, the vertex corrections do not sig- 
nificantly modify the bare-bubble results. 



representation, as in both methods have the same 

result. For impurity fraction z ranging from 0.5% to 1%, 
the importance of the vertex corrections is largely seen to 
be negligible, which implies that an analysis of the bare 
bubble results is sufficient. 

Figs. [51ITT] illustrate the vertex corrected thermal con- 



ductivities, 



in the main graphs, while the insets 



display the relative discrepancy with respect to the bare 
bubble thermal conductivities " ^bu — • Each is plotted 
as a function of the amplitude of the CDW, ip/i^c, where 
ipc indicates the maximal CDW for which the clean sys- 
tem remains gapless. We will postpone analysis of the 
character of the thermal conductivity until Sec. V C. 

To gauge the importance of the vertex corrections, we 
look first at Fig. [51 The inset indicates that the vertex 
corrections do not signifigantly modify the bare bubble 
thermal conductivity. Although their importance grows 
somewhat with increasing t/j, the correction is still slight. 

Next, Fig. [S] is used as a reference against which to 
consider the dependence of vertex corrections on scatter- 
ing potential, impurity fraction, and velocity anisotropy. 
The next three figures are the results of computations 
with each of these parameters modified in turn. By com- 
paring Fig. [9] with Fig. [8] we conclude that the vertex 
corrections become more important when the scattering 
potential is peaked in fc-space, but are unimportant for 
potentials that fall off slowly in /c-space. 

Fig. [5] and Fig. [TO] correspond roughly to the largest 
and smallest z for which these calculations are valid. 
Comparison of these two figures, as well as that of in- 
termediary values of z (not displayed) indicates that the 
relative importance of the vertex corrections is indepen- 
dent of z. Nor does increasing the velocity anisotropy 
affect their importance, as seen by making a comparison 
between Fig. [5| and Fig. [TT] 



B. Clean Limit Analysis 

It is of great interest to consider the behavior of the 
thermal conductivity in the clean (z 0) limit. Because 
the thermal conductivity is composed of integrals over 
p-space of functions which become increasingly peaked 
in this limit, there exists a sufhciently small z beyond 
which it is not possible to perform the requisite numer- 
ical integrations. However, it is still possible to obtain 
information about this regime. To that end, we will ex- 
amine the form of the bare-bubble thermal conductiv- 
ity, and consider the z ^ limit. As we shall see, this 
will enable us to determine the value of -0 at which the 
nodal approximation, and hence this calculation, is no 
longer valid. Additionally, a closed-form result for the 
thermal conductivity in the z — )■ limit is obtained for 
the isotropic {vf = va) case. The bare-bubble thermal 
conductivity, identical with setting A ^ in Eq. (|5T|) , is 
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kg Vf + Vl 
3 VfV2 ' 



d^p Ni + N2 



2tt 



D 

2\2\ 



ei = Cfc 



£2 = f-k+G 



77„,A (V' - B^Y{{ei + - (Ai - A2)") + ((Ai + A2)' - (ei - £2)') - 4Si(V' - B3)(eiAi + 62^2)) 



Ni = A{{A + B + el + AlY + {A + B + 4 + A 
N2 

D = \{A + B + ei + A^)(A + S + + A^) - B(^(ei + £2)" + (Ai - A2) 
+4Si(^Bi(eie2 - A1A2) + (V' - -B3)(eiA2 + e2Ai; 



Ai = Afc 

A2 = Afc+G 



(52) 



where A = Tl and i? = (ip - B^Y + Bf. Since the 
results of Section UlI Bl indicated that Tq ^ exp (— ^) and 
_Bi, i?3 '--^ z, in the z ^ hmit, A — > much faster than 
Bl ^ OT B3 ~^ 0. Therefore in taking the z — > Hmit 



we will first let ^ — > to obtain a result still expressed in 
terms of Bi and B^. The denominator can be rearranged 
as 



D = (^{A^ + A{2B + el + Aj + el + Aj) + where 

/ = B^ + iel + Al){el + Aj) - 2B{eie2 ~ A1A2) + 4(Bi(eie2 - A1A2) + - B3)(eiA2 + e2Ai)) 

= ((£162 - A1A2) - {2Bl -B)y + ((eiA2 + e2Ai) + 2Bi(^ - B3))' (53) 



We are thus considering, in the limit that A ^ 0, an 
integral of the form 

/ d^P ^-^^P^ 2 (54) 

J {Ah{p) + /(p))' 

Note that any nonzero contribution to this integral must 
come from a region in p-space in which /(p) — 0. We 
will consider separately the isotropic case {vf = va) and 
the anisotropic case {vf > va)- 



I 

1. Isotropic Case 



For the special case where Vf — va, it is possible to 
calculate the integral of Eq. ([52]) exactly, by taking the 
A ^ limit, and choosing another parametrization. The 
coordinates qi = Ck - efc+g and 172 = Cfc + efe+g - 1, 
have their origin located at the midpoint of the white 
and gray dots of Fig. [1] Using these coordinates, in the 
yl — > limit we find that the elements of Eq. (|5^ become 



iVi = 2Ai^B^ + B{q^ + 1) + -{q^ + If + - 

N2 = 2r/™A((^ - B3)'((ei + £2)' - (Ai - A2 f ) + BI{{Ai + A2)' - (ei - £2)') - 4Bi(^ - B3)(eiA2 + e2Ai; 

= 2r,„((V^ - i?3)'(g2' + 2^2 + 1 - ql) + Bliql - 2q2 + 1 - qf) - AB^i^J - Bs)i2ql - q^ - 1)) 

= 27jrnA[{2ql ~q^ + 1)((^ - Bsf - 2Bi(^ - B3) + Bfj + 2g2((V' - B3)' - Bf^ + ABi{-^ - B, 

D 



2A{i + B- 2Bi(^ - S3)) + (92 - (^' - Bff^ + _(q2 _ _ 4^^(^ _ 5^))^ 



(55) 



r 



Now the part of the denominator not proportional to A, In qi/q2 coordinates, these are the equations of a hori- 
thc /-term, is zero when 



92 = (^-53)2-5? and ^ 1 - 4Si(^ - S3). (56) 
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FIG. 12: Illustrated is a schematic view of the line and cir- 
cle whose intersection determines whether gapless excitations 
remain, for the isotropic case {vf = va)- The left figure indi- 
cates the situation in the absence of charge ordering, that is, 
for i/> = 0, where the radius of the circle is 1 and the line lies 
on the horizontal axis. The right figure indicates the situation 
at i/) = i/)*, when charge ordering is such that the excitation 
spectrum becomes gapped. In the clean case, the evolution 
corresponds to moving the line past the circle. With self- 
consistent disorder, the radius of the circle and height of the 
line are both functions of tp; in each instance, this construc- 
tion can be used to determine the value of ^ at which the 
quasiparticle spectrum becomes gapped. This value of i/> is 
referred to as tp* in this paper. 



zontal line and a circle, which must intersect for there to 
be a nonzero contribution to the integral, since each term 
is positive definite. In the simplified disorder treatment 
of Ref. 113 for which Bi = ~ and Tq = constant, 
these constraints simplify to q2 — tp"^ and — 1, so that 
no contribution occurs when ^ > 1 (Note that as in the 



numerical analysis, being an energy, is measured in 
units of ipc)- With the self-consistent treatment of disor- 
der, there will likewise be a sufficiently large value of 
beyond which the line and circle no longer intersect; we 
will call this value ip* (see Fig. [T2|) . We interpret ip* 
the point beyond which the system becomes effectively 
gapped. This is consistent with the exact result found by 
computing the eigenvalues of the completely clean hamil- 
tonian (as "0* = V'c m that case). 

In Sec. nil Bl it was determined that Bi ~ biip and 
B^ ~ &3'0j where bi and 63 depend on the remaining 
parameters of the model. Using this approximate form 
for Bl and B3, the condition for the maximum ip for 
which the constraints of Eq. are satisfied, 

I _ 4Bi{^ - B3) = (^{^ - Bsf - Biy , (57) 
indicates that 

±((1-&3)t6i)' 



(^{i-b3-bi){i-bs+bi)y 



(58) 



Since ip*'^ > 0, we find that for Vf = va, 



1 



(59) 



1 - 63 + bl 

We now proceed with the calculation of the clean-limit 
thermal conductivity. Substituting the conditions of 
Eq. ([SS]) into Eq. we find that the numerators be- 

come 



AA 1-2^1(^-53) l + B-2Bi(0-B3) 



N2 = AT^„Ml + B-2Bi{4,~B^)]{[{'4,~B:,f -Blf + 2Bi{4,-B- 



(60) 



both of which are independent of q, so that the clean The details of this integration are reported in Appendix 

limit result hinges upon the integral [Cl with the result 

/=/— ^ 2' (61) 

■I '^^ (fclA+(g2-fe)2 + i(q2-fc3)2) 

1= ^ . (63) 

where 2fci-\//c3 — fc^ 

fci = 2(1 + B - 2Bi{i' - B:i)) 

^2 — ~ B^) — Bl -y^g gg^jj j^Q.^ write the anisotropic clean limit thermal 

^3 = 1 - 4:Bi{tp - B3). (62) conductivity 
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1 - 2Bi^ + ?7„, ([(^ - B^f - BlY + 2Bi(^ - B3) , , 
^1 - 4Bi(V^ - B3) - [(t/- - ^3)2 - BlY ^ 

= ^1 - 4Si(7A - B3) - [(V' - ^3)2 - e(i-4i3i(V'-B3)-[(V'-S3)'-i3?]') 

- . 'M'^''".'''rr''p'.2 P.I. e(l-4S,(^-S3)-[(^-S3)^-S,T), (64) 

y/l - ABiiiJj - B3) - [iij - B-iY - BlY ^ ' 



I 

where the function is the Heaviside step function. Us- we are able to rewrite the dimensionless conductivity in 
ing the definition for '0* found in Eq. ([SQ]). and defining terms of parameters easily extrapolated from SCBA cal- 
culations 

I 




in which form it is clear that the thermal conductivity 
vanishes for ip > ip*. 



2. Anisotropic Case 

For the case of anisotropic nodes, Vf > va, the integral 
of Eq. ([5^ becomes intractable. However, it is still pos- 
sible to predict Using the same 91/(72 coordinates, 
the /-part of the denominator is again a sum of two pos- 
itive definite terms. Again, the only contributions to the 
clean-limit thermal conductivity arise when / = 0, which 
again provides two equations 



(67) 



where 



2/3 



1 



/34- 



-^l-{f3^-l){{^-Bsr-Bf) 



R = 



(l--(/3-l))2-4i?i(V'-i?3). (68) 

P 



This defines a hyperbola and a circle, again parametrized 
by ij}- One instance of this is depicted in Fig. [T31 The 




FIG. 13: For generally anisotropic Dirac quasiparticles, the 
construction used in Fig |12l is modified to contain a hyper- 
bola and circle. When these no longer intersect, the excita- 
tion spectrum becomes gapped. Illustrated is the construction 
for scattering parameter values {Vi, R2, R3} = {110, 0.9, 0.8}, 
impurity fraction z = 0.01, and with Vf = ivA- For these 
parameters it was determined that the value of at which 
the spectrum becomes gapped is given by tp* = 0.32ipc. 



value of at which these equations no longer have a 
solution is Tp* ■ The computed values for ip* are included 
for comparison in the plots of thermal conductivity in 
Fig. [H and Fig. [13 
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0.2 0.4 0.6 0.8 1 

FIG. 14: Effects of disorder on the charge-order-dependence 
of the bare- bubble thermal conductivity, isotropic case {vf = 
va)- Note how an increase in the impurity fraction, z, broad- 
ens out the peak in the conductivity. As the disorder becomes 
sufficiently small, the computed conductivity (triangles and 
squares) attains a limiting value that closely agrees with the 
closed-form clean-limit results of Eq. (|66p (shown with solid 
lines). The thermal conductivity obtained by simply letting 
E ^-iFo (as in Ref. |47j i) is shown with dashed lines. The 
effect of the self-consistent disorder is to renormalize the ef- 
fective tp at which the thermal conductivity vanishes (from 
4>c to ipc)- Here, we have considered short-ranged scatterers 
{Vi,R2,R3} = {110,0.9,0.8}. 



C. Effect of Self-Consistent Disorder 

Satisfied that vertex corrections are of little impor- 
tance, we set about analyzing the form of the thermal 
conductivity by studying the bare-bubble results. Ther- 
mal conductivity k was computed for /3 = \/vJJv^ val- 
ues of 1, 2, 3 and 4 (that is, for u//wa = 1, 4, 9 and 16). In 



0.6 0.8 

FIG. 15: Effects of disorder on the charge-order-dependence of 
the bare- bubble thermal conductivity, anisotropic case {vf = 
16va)- The effect of disorder is the same as in the isotropic 
case, which is to mix gapped and gapless states, smearing the 
peak in Hyy across the renormalized nodal transition point. 
It is interesting to note that for this anisotropic case, tp* 
is significantly smaller than tpc- Again, we have considered 
short-ranged scatterers {Vi,-R2,-R3} ~ {110,0.9,0.8}. 



Fig.[T3]is presented a representative plot of k for Vf = va- 
The clean limit prediction for k (Eq. ((66)) ) is computed by 
fitting bi and 63 from the self-energy calculations. These 
clean limit predictions are then plotted on the same graph 
with the numerical results of the thermal conductivity for 
the same parameters. In addition, the clean limit results 
of the simpler disorder model of Ref. [13 are also shown for 
the Vf — Va case. Increasing disorder broadens the peak 
in near ^p*. For z — 0.005, the numerical compu- 
tation is already almost exactly given by the clean limit 
results, while for z = 0.009, the features of the conductiv- 
ity are nearly totally smeared out, as seen in Fig.[T31 In 
this figure, the value of ip* given by Eq. (|59p is indicated 
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with an arrow. 

For Vf > v^, the thermal conductivity has the same 
characteristics as for Vf = ua, except that ip* is gen- 
erally smaller for larger /3. The numerically computed 
thermal conductivities for the case of /3 = 4 are shown 
in Fig. 1151 In this figure, the value of ip* is computed 
by determining the largest value of -0 for which Eqs. (j67p 
have a solution, and is indicated with an arrow. It is clear 
from these graphs that the self-consistent disorder renor- 
malizes the amplitude of charge density wave at which 
the thermal conductivity vanishes, and that the amount 
of renormalization is heavily dependent on the velocity 
anisotropy ratio, and varies only slightly with changing 
impurity fraction. 

VI. CONCLUSIONS 

The work described in this paper investigates the low 
temperature thermal conductivity of a d-wave supercon- 
ductor with coexisting charge order in the presence of 
impurity scattering. We improve upon the model stud- 
ied in Ref. 113 by incorporating the effect of vertex cor- 
rections, and by including disorder in a self-consistent 
manner. Inclusion of vertex corrections does not signif- 
icantly modify the bare-bubble results for short range 
scattering potentials. The role vertex corrections play 
increases somewhat for longer range scattering poten- 
tials, in particular as the amplitude of charge ordering 
increases. Nonetheless, for reasonable parameter values, 
the inclusion of vertex corrections is not found to signifi- 
cantly modify the bare-bubble results. This opens up the 
possibility of doing bare-bubble calculations for models 
with different types of ordering. 

Our analysis determined that for self-consistency, it 
is necessary to include off-diagonal (in extended-Nambu 



space) terms in the self-energy. As the charge ordering 
increases, the off-diagonal components become more im- 
portant, and are found to dominate the self-energy in the 
clean limit. We also find that the zero-temperature ther- 
mal conductivity is no longer universal, as it depends on 
both disorder and charge order, rather than being solely 
determined by the anisotropy of the nodal energy spec- 
trum. 

In addition, inclusion of disorder within the self- 
consistent Born approximation renormalizes, generally 
to smaller values, the critical value of charge ordering 
strength ip at which the system becomes becomes effec- 
tively gapped. This renormalization is seen in the calcu- 
lated thermal conductivity curves, and depends primarily 
on the impurity fraction z and velocity anisotropy w/ / va ■ 
For larger Vf/vA, the renormalization can be significant, 
which may indicate that the calculated effects could be 
seen in low-temperature thermal transport even in sys- 
tems with relatively weak charge order. 
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APPENDIX A: SELF-CONSISTENT GREEN'S 
FUNCTIONS 

Here are the Green's functions that fulfill the self- 
consistent Born approximation. The superscript ^'^^ refers 
to the fact that 3 successive applications of our self- 
energy scheme were necessary for self-consistency, as is 
explained in Section III. 



+4 (^flUA + /3pi)(^c + /3P2) - ^PiP2)) - fs^HA + f3pi)pi + (V'c + f3p2)P2)^ 

-if 2 + /3')((2^c + PiPl +P2)f + -^iP2 -Pi)') 
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41(^ 



Pl,P2) 
Pl,P2) 
Pl,P2) 
Pl,P2) 
Pl,P2) 
Pl,P2) 

Pl,P2) 

Pl,P2) 
Pl,P2) 
Pl,P2) 
Pl,P2) 
Pl,P2) 
Pl,P2) 
Pl,P2) 



= -h -ft + i^c + f3p2Y + -^Pi + n + /a 



/32 

-pP2^-fl + {i^c.+PP2f + {-^Plf^ '--^ 



^Pl(/|-/l)+2(V'c + /?P2)/2/3 



-(^, + /3pi) ( + (V-c + Pp2f + {^Plf^ + (V'c + /?P2)(/3' - /I) + |pi/2/3 



A ( /alZV-c + /3bl +P2)) + /2^(P1 +P2) 



/2 (^/l' - (V^c + /3pi)(V'c + /3P2) + ^PlP2 - /I - /l) + h ((^c + + (V'c + /3P2)P2) 



/3 (^/l - /2 - /a + (^c + Ppi){4'c + /?P2) - -^PlP2 

41(^;pi.P2) 

41(^;pi'P2) 

-41(^;pi'P2) 

41('^;pi'P2) 

41('^;p2,Pi) 

41(^;p2,pi) 



/2 ((V-c + f3pi)pi + (V'c + PP2)P2) 



(Al) 



To obtain the retarded Green's function G'^°*(ijj) from 
the above we set 



/2 = Siricc) 

/3 = V' + Sg°3*M 



(A2) 



For the retarded Green's function G'^°*(w + il), we set 
LD — !■ w + and for the advanced Green's function 
G^^^{lo) we set ER<=' ^ E^'*^ by taking the complex 
conjugate. 



APPENDIX B: CUTOFF-DEPENDENCE OF 
SELF-ENERGY 

Here we note that the self-consistent Born approxima- 
tion, when applied to the nodal Green's functions used 
in this paper, produces a self-energy that contains a log- 
arithmic divergence, and therefore has a prefactor that is 
proportional to the momentum cutoff, set by the size of 
the Brillouin zone. By contrast, the thermal conductivity 
has no such dependence, and is therefore a truly nodal 
property. One difficulty this introduces is that the pref- 
actor of the self-energy is sensitive to our choice of coor- 
dinates. As the location of the nodes evolves with charge 
density wave order parameter tp, computations are nec- 
essarily performed in a different local coordinate system 



(than one centered about a node itself). This coordinate 
shift in the pi direction introduces a constant Tia^ term, 
even in the il) — Q instance (whereas using node-centered 
coordinates, the anti-symmetric integral is found to van- 
ish). In the V' = case, a shift of e corresponds to the 
integral 



/ = 



dpi 



Pa+t 



Po 

dp2 — 

-pa Pi 



Pi 



pI 



^ 



(Bl) 



The result is 7re, which matches the discrepancy. We 
therefore subtract off the = value of E^a; the re- 
sults shown in Fig. [3] reflect this recalibration, as do the 
subsequent iterations of the self-energy calculation. 



APPENDIX C: CALCULATION OF CLEAN 
LIMIT INTEGRAL 

For the clean limit of the thermal conductivity we need 
the integral 



d2g 



A 



(CI) 



in the limit A ^ 0. With the substitution 

Xi = xcosd — qi — 1 X2 = xsvo-O = q2 (C2) 
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the quantity Yi = {q2 — + \{<1^ ~ ^i)^ becomes To simplify the angular integrand, we get rid of the sin0 

term by shifting d ^ + a. Then, the last two terms of 



2 , /l-fc3\2 ,1-^3 2 , 2 Eq.Elbecome 



+ (.t2 + 1 - fc3)a;cos6'- 2a;A:2sin6'. (C3) 
I 



008(6/ + a) — A;2 sm(6^ + a) = ( cos a — K2sma)cost/ 

— ( sina + fc2 cosa) sin0. (C4) 



We set the coefficient of the second term on the RHS of Eq. (|C4p to 0, so that the first term becomes 

" k^V 2 + 'J smacos^ = _(^_ + -^x^ + i^^) + ^Ij cos0, (C5) 

where the RHS of Eq. (jCSp is obtained by setting sin a = k2/r, where r — r{x) is an undetermined function of x. 
With this substitution, Eq. (jC3[) becomes 



X'^ 1 — ^3 9 /I — ^3x9 ,9 9 2xfx'^ 1 — ^3 9 /I^fc3x9 , 9\ //, N 

( (^^±1^)2 + kf] { 1 + f cos(0 + a)) 

2 

^(l + a^ -2acos(6l + a)), (C6) 



4 2 

, X"^ 1 /l3 o , 9\ 2X 



where 



we find that 



a;2 + 1 - K \ 2 



fc| and 



a=-. (C7) 
r 



Then, defining 7 = k\a^ jx^^ the integral of Eq. (|C1 
becomes 



d2x A 

^/ciA +f^(l + a2- 20008(6* + 



27r a;4 



A d6l 



° (^^7+ 1 + a2 - 2acos(6' + Q;) 



^8) 



after shifting 6* — > 6* — a, and noting the evenness of the 
Q integral. The integral is found in standard integration 
tables^, and noting that (1 ± of' + ^7 > 0, we obtain 



/ = 



da; A-k(\ + a^) 



;ii±P((i-«,'..,)--,(c.) 



/o 27rx3 (1 + a) 
Since in the limit that ^ — > 0, 

A 2 



((1 + 0)2 + ^7) 



3/2 



7 



(5(1 -a), (CIO) 



/ = 



da; 



4fci 



(5(a- 1). 



(Cll) 



1,2 

«,2 



Making the further substitution y = (a;2 + 1 — k^)/2^ 



00 

1 -^3 

2 

00 

1 -^3 

2 



dy 



1 



2fci j/2 
dy 



A;| 



2;/-(i-fc3) 

(y2 + fc2)2 



- 1 



4fci 



A:2 - j;2 + _ 



:(s{y-y+)+S{y-y^)j, 



where 



y± 



l±\ k 



(C12) 



(C13) 



are the intersections of the curves j/2 fc2 g^^^^^ 2y — (1 — 
fca). It is easily verified that both y^ and y_ are in the 



range of integration 



l-k3 



, 00) (y_ just catching the lower 



bound when = 0). Then expanding the denominator 
of Eq. (|UT2|1 using Eq. (|UT3l) . we find 



18 



H y± + y±ii - h) 



2x1 h - kl 



i + fc. 



(C14) 



so that 



/ = 



1 1 /l + fe3 + 2A/fc3-fc^ 1 + fcg - 2^fc3 - k. 



2fci 2 v^fca - fc2 V 1 + + 2^fc3 - ^ l + h- 2^fc3 - ^1 
1 

2fci\/fc3-fc^ 



(C15) 
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